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Feynman Path Integration in Optics: Numerical Example * 


1. Introduction. 


In the following article we want to show how to solve the problem of propagation of 
light in a planar waveguide structure with a graded refractive index. The method here is a 
consequence of Feynman path integration formalism [1]. This formalism was developed for 
quantum ,mechanics and can be applied to solve the typical quantum mechanical 


problems [1,2]. 


#) Work supported from RR102 Fund. 
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2. Schrodinger equation in paraxial optics. WA) 

In the case of light propagation in media with graded index profile the Schrédinger— 
ike equation appears in the so called paraxial optics [3,4]. In order to have such a 
situation consider the planar waveguide directed along ғ axis. The coordinate x is 
perpendicular to the propagation direction and measures the lateral dimension of the 


waveguide. The Schrédinger equation for this case reads [4]: 
3 
i 2- Æx s) a [= k ү: + а(х) E(x,2) = H Е(х,з). (1) 


Here n,(x) is the deviation of the refractive index n(x) from its constant value ng: 
n(x)=np+n,(x), is the wave lenght of propagated light and Е(х,2) is the electric field in 
the waveguide. In this equation the role of the Planck constant is played by the reduced 
wavelength A/27 of propagated light. 
To solve this equation one can use the Green function G(x,x';z) which is defined as 
a propagator of the point signal as follows: - 
[н - 52 G(x,x!;2~Z9) = “ih Ha) &z—2'). (2) 
The Green function G(x,x!;z) is related to propagation constants {8,2} through the relation 
known from the analog situation in quantum mechanics. The formula reads: 
С(х,х!;2) = bite u mê) ехр(18 12) а(х), (3) 
where u m(*)=<x|m> is a normalized eigenfunction and Pm is the eigenvalue of the Н 
operator defined in eq.(1): 
H|m>=—-{A/2z)6_|m>. (4) 
By the use of the orthonormality of the modal eigenfunctions um), one can find the 
useful relation between С and {£p}, which does not exploit the explicit form of these 


functions. One has 


ІН) This chapter and Chapter 3 follow closely the material of paper [I] by Hawkins and 
serve as а guide for understanding Chapter 4. 
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ш 
о 
Tr G(x,x';z) = Sf ох») ах = Ў expli з). (5) 
—@ m=0 
In the next step we use the analitical continuation of the Green function G(x,x';z) 
by the substitution 7 = —iz. If 7 is large and positive then the largest contribution to the 
sum in equation (5) comes from the lowest eigenvalue Af, of Н operator. Thus one has the 
relation 


By =—-lim 1 ln Tr G(x,x';i7), (6) 


T> wo 


which is the known Feynman—Kac formula [2]. In the same limit of large т one can also 


obtain the intensity distribution function 


ана О(хх дт 
fos) = lim т (т) 


3. Тһе evaluation of the Green function. 

In order to evaluate the Green functions needed for calculation of the propagation 
constant Ay of the fundamental mode and the intensity distribution function 15 опе сап use 
the short time approsimation [3,5—7] and exploit the group property of the evolution 
operator exp(—tz/). In the case of light propagation equation the short time corresponds 
to short distance. For the short distance slice 2е one has 

@ 
G(x,x'2ie) = f Сб(хх"дє) GO" xie) dx! (8) 
“ 


-Ф 


If we chose the final distance equal to 7 then the successive iterations of this equation give 
the Green function corresponding to 7. In the following we take =N e where 2N is the 
convenient number of distance slices in the distance interval т, each equal to e. Thus, we 
have made the distance discretized for our numerical purposes. If N ( or 2% is large then 
this approximation converges to the "smooth" Green function G(.,.;i7). 


The short distance propagator G(x,x';ie) is given by [3,5,6] 
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G(x,x';i¢) = aje A) |а 
s exp{— 2° ахуа е) expl- ya). (9) 
The integral іп eq.1 сап be calculated by using rectangular quadrature approximation of 
the form 
M 
С(х х2) Е ч 2) С(хұх.3е) Сх худе), (10) 
where x,€(-C,C) and j=1,....M. The parameter A із the griding distance of x coordinate 
and is small enough to make the quadrature approximation as good as possible. The cut off 
parameter C is finite in calculations. However its value is larger then the lateral measure of 
the considered waveguide. As it was pointed out іп [3] the convergence of the iterating 
procedure is obtained when both parameters А and є fulfill the condition ny A*/re < 0.1. 
This last criterion was also tested in our calculations. 

The iteration procedure described above can be make very fast if one exploits the 
properties (reality and symmetry) of the Green matrices G(ie) for the infinitesimal distance 
є and G(i7) for the finite distance 7. From equation (10) one has 

Gin = a- [GG9) 2“ (1) 
Diagonalizing геа] and symmetric matrix С by a proper matrix S one finds 

р(їт) = 57106(17)8 = sta ае sea? ID? Gg), | 12) 

where D(...) is the diagonal matrix with the eigenvalues of G(...). These results together 


with the approximation (10) give for the trace operation (eq.5) 
APTON 
Tr G(ir) = A Edik Ў (13) 


where ау аге the diagonal elements of D(ie). 


3. Example 1: Single strip. 
Consider the simple case of a strip waveguide with the graded refractive index in 


which the index deviation is 
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(х= саана (14) 

1 cosh“ (2х/һҺ) 
where бп їз the maximum deviation from ng and h is the lateral measure of the strip. The 
solution of the analogous quantum mechanical problem is discussed in Landau and Lifshitz 
textbook on quantum mechanics [8] and the expression for propagation constants {8} is 


given in the Kogelnik's paper [9]. The formula reads 


‚ 2 
fa =~ AR (s—m)*, (15) 
where 
2s={[1-2 (24822. іу, (16) 


and ш is the mode number. We will use the following values of the waveguide parameters: 


ng=2.2, бп--0.006, corresponding to the Ti diffused in LiNbO, [3,9]. 


Table 1. Propagation constant for zero mode 8 (А-1.3 ша). 


k (pm | Ay(calc) Aplexact) 
6 —0.01903 —0.01903 
Т - -0.02020 —0.02020 
8 0.02112 —0.02112 


ж) This result coincide with that in Table 1. of paper |4). 


Figure 1 shows the index profile and the calculated intensity distribution Ip(x) for 
the fundamental mode in the case of strip with lateral dimension h=6yum and the typical 
length of the propagating light, A=1.3 шп. Exact and calculated values of By for different 
widths h are shown in Table 1. Notice a very good agreement of the calculated values of 
propagation constants with the exact ones. This example serves also as the test of our 


numerical procedure. 
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5. Example 2: Double strip. 

As a more complicated example of application of the Feynman path integration let 
us consider the two strip planar waveguide with completely symmetric index profile. The 
index profile is modelled by smoothly joint parabolic curves or by two Gauss functions. 
The parameters listed in Figure 2 as well as in other figures in this paper are chosen so as 
to be close to realistic situation like this in Chapter 3. The results of numerical calculations 
are shown in Figures 3 — 10. By the use of different models and various parameters one can 
show the basic properties of the solutions of Schrödinger equation (1). 

Figures 2 to 6 show the behaviour of the intensity function for various widths and 
heights of the index profile and for different distance between strips. For the symmetrical 
index the intensity distribution is also symmetrical. In the case of well separated strips the 
disturbed symmetry causes the signal to disappear either in the lower strip or in the 
narrower one. For more overlapped strips the intensity distribution is located in one of the 
strips and shows a hump at the side of the other strip. 

The similar situation can be observed for the refractive index built of two Gauss 


functions. The deviation of the refractive index is 


sale (xx) én (х-х2)7 
= a xp[- 20} 1+ Boat 54 1. (18) 


Неге the standard convention for parameters із used. The results are shown іп Figures 7 to 
10. The interesting results for this case are those for asymmetrical widths and/or 
asymmetrical heights of the refractive index profile. The propagation constants for . 


fundamental mode are given in each figure. 


6. Discussion. 
The above examples show that the method of path integration is effective in the one 
dimensional case and for various kinds of refractive indices. Its application is independent 


on the symmetry property of the index profile. The similar procedure can be applied to 
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determine the propagation constants for higher order modes. 

The path integration can also be realized by using the Monte Carlo method [6]. 
However the calculation times in this case are one order of magnitude larger as compared 
to these of the present method. 

It seems that the described procedure can be succesfully applied to two dimensional 
cases. However the large matrices entering the calculations will considerably slow down the 
effectivness of the matrix multiplication algorithm and the time of calculations will be very 


large. 
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Figure 6. 4 
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Figure 7. 
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Figure 9. 
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